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I.  INTRODUCTION 


The  purpose  of  this  report  is  to  describe  a six  degree  of 
freedom  missile  trajectory  program  in  sufficient  detail  for  the  user  to 
be  able  to  run  the  program  and/or  make  modifications  to  fit  his  specific 
needs.  This  program  utilizes  a nonrolling  missile  coordinate  system 
to  speed  up  run  time  when  high  (greater  than  15  cps)  roll  rates  occur. 
Several  options  are  available  such  as  directional  control  force,  fold- 
out  fins,  and  print-out  of  steady-state  angles  and  velocity  angular 
errors  at  boost  burn  out. 


The  program  is  in  modular  form.  All  variables  to  be  integrated 
are  stored  in  an  array.  A subroutine,  DESUB,  consists  of  the  twelve 
simultaneous  differential  equations  describing  the  behavior  of  these 
variables.  A general  RUNGE-KUTTA  subroutine  is  utilized  to  integrate 
this  set  of  equations.  This  report  includes  a listing  of  these  equa- 
tions, programdescription  and  flow  charts,  definitions  of  all  symbols 
used,  and  a complete  program  listing. 


II.  programdescription 


Equations  of  Motion 


This  program  was  written  to  simulate  the  following 
equations  of  motion  to  describe  a free  rocket  trajectory  from  ignition 
to  impact. 


(T  - 1/ 2pSC  u )/m  + rv  - qw  - gSIN  (QE  + 9) 


ru+  (l/2pSu  C 


gCOS  (QE  + 9)  + qu  +(-l/2pSU  C, 
+ T8  SIN  («  + « ) - K®)/m 


r - (l/2pSdu  1 C. 


l/4pSd  ruC 


(l/2pSdu  lgC 
+ TS  dl  SIN 


l/4pSd  quC 


p a (1/ 2pSdu 


S - l/4pSd  puC  + l/2pSdu  Cin)/I  + P 


x - uCOS  (QE  + 0)  + wSIN  (QE  + 0)  COS*  - vSIN\(r 


y - uCOS  (QE  + 0)  SINi|r  + wSIN  (QE  + 0)  SINy  + vCOS* 

z - wCOS  (QE  + 0)  - uSIN  (QE  + 0) 

• 

0 - q 

t - r/COS  (QE  + 0) 

• 

♦ - p 

If  a constant  roll  rate  or  a programmed  roll  rate  is  to  be  imple- 
mented, this  can  be  done  by  setting  p « 0 and  £ “ ^gpj.n*  ^h*-8  *8  noted 

in  the  program  listing.  In  post  flight  analysis,  frequently  roll  history 
can  be  acquired  from  telemetry  data.  This  can  be  stored  in  tabular  data 
in  A ^ and  utilized  in  the  preceding  manner. 

These  equations  were  converted  to  Fortran  code  by  using  the  Fortran 
names  listed  in  Table  1.  These  equations  are  calculated  in  the  follow- 
ing coordinate  systems: 


Missile  coordinate  system. 


ft- 


B.  Program  Outline 

Solving  the  equations  of  motion  requires  current  values 
of  necessary  atmospheric,  aerodynamic,  and  missile  physical  character- 
istic data.  The  program  is  set  up  to  read  these  data  in  by  way  of 
NAMELIST,  This  is  one  of  the  easiest  methods  of  inputting  data.  An 
example  of  a data  card  for  a constant  is 

x = 75., 

This  is  done  with  no  regard  for  format  except  that  only  columns 
2 through  72  may  be  used.  If  x is  an  array  of  five  numbers,  it  is 
read  as: 

x (1)  - 6.,  4.,  1.,  17.,  0., 

This  can  be  continued  onto  the  next  card  is  necessary.  All  Fortran 
names  used  in  the  program  are  defined  in  Table  2.  Program  inputs  are 
listed  in  Table  3. 

Next,  all  variables  to  be  integrated  as  a funtion  of  time  are 
stored  in  an  array,  VARTB.  VARTB  (1)  is  time.  The  remainder  of  the 
array  is  ordered' consistent  with  the  order  of  the  corresponding  differ- 
ential equations  in  DESUB. 

The  current  values  of  atmospheric  data  are  now  calculated  using 
the  following  equations: 

G = Gq  [l/(l  - z/ 20855531. 5)]  2 
A * -z/(l  - z/20855531.5) 

If  A < 36098.2388, 

T - Tq  + (-0. 00356616  A) 

P =■  PQ  (t/t0)0,1633654Go 

p - 28.966P/49719.6257T 

Vs  -^2403. 075  T 
If  A > 36089.2388, 

T - 389.9878 

_ ,,,  -0.00000149G.  (A-36089.2388) 

P ■ 472.6765  e 0 
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P - 28.966  P/ 19390047. 44 

V - 968.075 
s 

where  G - gravity 

T ■ temperature 
P “ atmospheric  pressure 
p = atmospheric  density 
V ■ speed  of  sound. 

The  subscript  0 denotes  the  value  of  each  of  these  variables  at  sea 
level. 

Next,  the  wind  components  are  transformed  from  the  earth  to 
missile  coordinate  system.  These  values  are  used  to  calculate  missile 
total  velocity  and  angles  of  attack. 

Current  values  of  missile  physical  characteristics  and  aerodynamic 
data  are  now  calculated.  The  subroutine  XINT  Is  utilized  for  this 
purpose.  This  Is  simply  a linear  interpolation  routine  with  one  inde- 
pendent variable.  Logic  is  employed  in  this  section  to  determine  if 
the  missile  is  in  boost  or  coast  phase  and  whether  fins  have  been 
unfolded. 

Now  the  fourth-order  RUNGE-KUTTA  subroutine,  SOLVDE / S S IMDE * is  called. 
The  integration  time  step,  DELT,  is  an  input  to  this  subroutine.  A rule 
to  follow  in  determining  the  Bize  of  the  time  step  is  not  to  allow  any 
angle  to  change  more  than  7°  to  8*  during  one  step.  This  does  not 
include  the  roll  angle  unless  a thrust  misalignment  case  is  being  run. 

The  variable  DELMIN  is  the  minimum  the  program  will  allow  DELT  to  be 
changed  to.  This  is  done  in  SOLVDE/SSIMDE.  DELMIN  can  be  set  to  the 
same  size  as  DELT.  Other  inputs  to  this  routine  are  VARTB,  the  table 
of  variables  to  be  integrated;  WWW,  the  dummy  working  array;  IC;  and 
N,  the  number  of  variables  to  be  integrated.  DESUB,  the  routine  con- 
taining all  the  differential  equations  to  be  solved,  is  also  required 
as  an  input  argument.  SOLVDE/SSIMDE  is  a two-purpose  routine  for 
integration.  It  has  the  capability  of  solving  one  nth  order  differ- 
ential equation  or  n first  order  simultaneous  differential  equations. 

After  returning  from  the  integration  subroutine,  several  variables 
are  calculated  for  output  and  printed  if  it  is  time  for  printing.  If 
the  missile  altitude  has  become  positive  (below  sea  level)  or  if  maximum 


Sellers,  William  R. , SOLVDE/SSIMDE,  Directorate  for  Management 
Information  Systems,  US  Army  Missile  Comnand. 
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simulation  time  has  been  reached,  the  program  searches  for  the  next 
case  of  data.  Otherwise,  it  loops  back  to  calculate  atmospheric  data 
and  hence  back  through  the  program.  Program  output  is  defined  in 
Table  4. 


C.  Subroutine  DESUB 

This  subroutine  contains  all  differential  equations  to  be 
solved  in  the  order  that  the  variables  are  stored  in  VARTB  in  the  main 
program.  The  subroutine  inputs  are  V which  corresponds  to  VARTB  in  MAIN, 
and  J,  the  number  of  differential  equations.  A computed  GO  TO  as  a 
function  of  J controls  the  flow  through  the  routine.  For  example,  if 
J = 1»  the  first  differential  equation  is  solved  and  control  is  returned 
to  the  integration  routine.  The  output  is  F which  is  the  calculated 
value  of  the  differential  equation  for  the  current  value  of  J.  Following 
is  an  example  of^a  simple  DESUB  for  solving  three  simultaneous  differen- 
tial equations,  V(l)  = x,  V(2)  = y,  and  V(3)  - z. 

SUBROUTINE  DESUB  (V,  F,  J) 

COMMON  X ,Y,Z 
DIMENSION  V(4) 

GO  TO  (10,  20,  30),  J 
10  F = X 
RETURN 
20  F = Y 
RETURN 
30  F = 7. 

RETURN 

END 


Fortran  Name 


Symbol 


VARTB  (1) 
VARTB  (2) 
VARTB  (3) 
VARTB  (4) 
VARTB  (5) 
VARTB  (6) 
VARTP  (7) 
VARTB  (8) 
VARTB  (9) 
VARTB  (10) 
VARTB  (11) 
VARTB  (12) 
VARTB  (13) 


Fortran  Name 

Symbol 

TL 

1T 

CLA 

Lc 

ROLLI 

I 

X 

XI 

Xy 

CLD 

ci& 

DCANT 

5 

c 

CLP 

°1P 

CLO 

C10 

PSPIN 

Pspin 

TABLE  2.  SYMBOLS 


Fortran  Name 

Definition 

TEMPO 

Temperature  at  sea  level 

PRSSO 

Atmospheric  pressure  at  sea  level 

GRAVO 

Gravity  at  sea  level 

THETA 

Flight  path  angle,  pitch 

IC 

Initial  condition  for  RUNGE-KUTTA  routine 

VARTB 

Table  of  variables  to  be  integrated 

TT 

Current  value  of  time 

WO 

Initial  velocity  component  of  missile 
in  pitch  plane 

U 

Initial  velocity  component  of  missile 
along  the  longitudinal  axis 

VO 

Initial  velocity  component  of  missile 
in  yaw  plane 

RO 

Initial  yaw  rate  of  missile 

QO 

Initial  pitch  rate  of  missile 

PROLL 

Initial  roll  rate  of  missile 

XR 

Initial  range  of  missile  in  earth 
coordinate  system 

Y 

Initial  lateral  displacement 
in  earth  coordinate  system 

Z 

Initial  vertical  displacement 
in  earth  coordinate  system 

THETO 

Initial  pitch  angle 

PSIO 

Initial  yaw  angle 

PHIIN 

Initial  roll  angle 

WTIN 

Initial  missile  total  weight 

WT 

Current  missile  weight 

XXTIN 

Time  of  first  print  out 

XXT 

Time  of  next  print  out 

KELL 

Index  indicating  whether  missile 
is  on  launcher 

GRAV 

Current  value  of  gravity 

TABLE  2,  (Continued) 
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Fortran  Name 


Definition 


TEMP 

PRSS 

P 

OM 

WY 

WY1 

WX 

WX1 

WU 

WW 

VTOT 

ALP 

ALPHA 

UM 

ALTFT 

XT 

T 

A 

PSPIN 

TTA 


Current  value  temperature 

Current  value  of  atmospheric  pressure 

Current  value  of  atmospheric  density 

Current  value  of  speed  of  sound 

Boost  phase  lateral  wind  component 

Coast  phase  lateral  wind  component 

Boost  phase  head  wind  component 

Coast  phase  head  wind  component 

Component  of  head  wind  along 
longitudinal  missile  axis 

Component  of  head  wind  along  vertical 
axis  in  missile  coordinate  system 

Total  missile  velocity 

Angle  of  attack  in  vertical  plane 

Total  angle  of  attack 

Current  Mach  number 

Current  altitude  (positive  up) 

Current  missile  thrust 

Time  table,  independent  variable  for 
thrust  table 

Thrust  table,  f(T) 

Roll  acceleration  due  to  spin  motor 

Time  table,  independent  variable  for 


AR 

TTA 


PSPIN  table 

PSPIN  table,  f (TTA) 

Time  table,  independent  variable  for 
ASPIN  table  (Note  both  ASPIN  and 
PSPIN  will  not  be  used  in  same 
simulation) 


AK 

ASPIN 

CLP 


ASPIN  table,  f (TTA) 

Constant  or  programmed  roll  rate 
Roll  damping  coefficient 
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TABLE  2.  (Continued) 


Fortran  Name 

Definition 

UG 

Mach  number  table,  independent 
variable  for  CLP  table 

AO 

CLP  table,  f(UG) 

CLD 

Roll  moment  coefficient  due  to 
fin  deflection 

UH 

Mach  number  table,  independent 
variable  for  CLD  table 

AP 

CLD  table,  f(UH) 

XK 

Control  force  for  system  with 
direction  control 

TTB 

Time  table,  independent 
variable  for  XK  table 

AL 

XK  table,  f (TTB) 

CMQ 

Pitch  damping  coefficient 

UE 

Mach  number  table,  independent 
variable  for  CMQ  table 

AN 

CMQ  table,  f(UE) 

CNA 

Normal  force  coefficient 

UB 

Mach  number  table,  independent 
variable  for  CNA  table 

AG 

CNA  table,  f (UB) , power  on 

CLO 

Roll  moment  coefficient  due  to 
wraparound  fin  effect 

UI 

Mach  number  table,  independent 
variable  for  CLO  table 

AS 

CLO  table,  f(UI) 

RNLM 

Reynolds  number  per  body  length 
times  Mach  number 

ZAL 

Altitude  table,  independent 
variable  for  RNLM  table 

AU 

RNLM  table,  f (ZAL) 

BL 

Missile  body  length  (ft) 

BLOD 

Missile  length  divided  by  diameter 
d/d) 
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TABLE  2.  (Continued) 


Fortran  Name 

Definition 

RN 

Reynolds  number 

CDF 

Coefficient  of  drag  due  to  skin 
friction 

TB 

Boost  burn  time 

SPIMP 

Specific  impulse 

TIMP 

Total  impulse 

TPROP 

Total  propellant  weight 

DELT 

Time  integration  step 

XII 

Initial  transverse  moment  of  inertia 

XIB 

Final  transverse  moment  of  inertia 

TTO 

Initial  time  for  table  values 

XI 

Current  value  of  transverse  moment 
of  inertia 

XL  I 

Initial  thrust  lever  arm 

TLF 

Final  thrust  lever  arm 

TL 

Current  value  of  thrust  lever  arm 

CLAI 

Initial  value  of  control  lever  arm 

CLA 

Current  value  of  control  lever  arm 

RDLLI 

Current  value  of  roll  moment  of  inertia 

ROLLA 

Initial  value  of  roll  moment  of  inertia 

ROLLB 

Final  value  of  roll  moment  of  inertia 

SL 

Static  margin 

UA 

Mach  number  table,  independent 
variable  for  SL  table 

AI 

SL  table,  power  on,  f(UA) 

CD 

Drag  coefficient 

UC 

Mach  number  table  independent 
variable  for  CD  table,  power  on 

AH 

CD  table,  power  on,  f(UC) 

UA 

Mach  number  table,  independent 
variable  for  SL  table,  power  off 

AJ 

SL  table,  power  off,  f(UA) 

TABLE  2.  (Continued) 


Fortran  Name 

Definition 

UD 

Mach  number  table,  independent 
variable  for  CD  table,  power  off 

AM 

CD  table,  power  off,  f(UD) 

UJ 

Mach  number  table,  independent 
variable  for  CNA  table,  power  off 

AT 

CNA  table,  power  off,  f(UJ) 

XM 

Missile  mass 

TFLD 

Time  of  fin  folding 

DELMIN 

Minimum  time  integration  step 

AX 

Velocity  along  the  x-axis  in  the 
earth  coordinate  system 

AY 

Velocity  along  the  y-axis  in  the 
earth  coordinate  system 

AZ 

Velocity  along  the  z-axis  in  the 
earth  coordinate  system 

DRAG 

Total  missile  drag 

ACC 

Total  missile  acceleration 

THETS 

Steady-state  pitch  angle  after  burnout 

PSITS 

Steady-state  yaw  angle  after  burn  out 

I BO 

Index  indicating  whether  to  print 
out  angular  errors  after  bum 

ZERR 

Pitch  angular  error  at  bum  out 

YERR 

Yaw  angular  error  at  bum  out 

XX 

Number  of  time-steps  per  print  out 

TMAX 

Maximum  time  - simulation  stops  after 
this  time 

LCS 

Index  indicating  last  case  of  data 

QE 

Quadrant  elevation 

S 

Missile  reference  area 

DIM 

Thrust  misalignment  angle 

PHIO 

Initial  thrust  misalignment  direction 

DX 

Missile  diameter 

I 


TABLE  2.  (Continued i 


Fortran  Name 

Definition 

ELL 

Effective  auncher  length 

DCANT 

Fin  cant  s :gle 

PCD 

Percent  elfective  drag  coefficient- 
utilized  in  error  budgets 

J 

Variable  number  to  be  integrated 

V 

Table  of  variables  to  be  integrated- 
peculiar  to  subroutines  DESUB  and  SOLVDE 

w,www 

Working  arrays  for  integration  routine 

N 

Number  of  differential  equations 

XX 

Array  of  independent  variables  for 
interpolation  purposes 

YY 

Array  of  dependent  variables  for 
interpolation  purposes 

X 

Current  value  of  independent  variable 

Y 

Current  value  of  dependent  variable 

TIME 

Current  value  of  time 

U 

Component  of  missile  velocity  along 
longitudinal  axis 

V 

Component  of  missile  velocity  in 
yaw  plane 

W 

Component  of  missile  velocity  in 
pitch  plane 

THET 

Pitch  angle 

PSIO 

Yaw  angle 

PHI 

Roll  angle 

Q 

Missile  pitch  rate 

R 

Missile  yaw  rate 

ACC 

Missile  acceleration 

VZ 

Velocity  along  the  z-axis  in  the 
earth  coordinate  system 

VY 

Velocity  along  the  y-axis  in  the 
earth  coordinate  system 

CNAI 

Constant  value  of  CNA  before  fins 
are  unfolded 

J 
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TABLE  2.  (Concluded) 


Fortran  Name 

Definition 

SLI 

Constant  value  of  SL  before  fins 

are  unfolded 

VARTB  (1) 

Current  value  of  time 

VARTB  (2) 

Current  value  of  W missile  velocity 

component 

VARTB  (3) 

Current  value  of  U missile  velocity 

component 

VARTB  (4) 

Current  value  of  V missile  velocity 

component 

VARTB  (5) 

Current  value  of  R,  missile  yaw  rate 

VARTB  (6) 

Current  value  of  Q,  missile  pitch  rate 

VARTB  (7) 

Current  value  of  missile  range 

VARTB  (8) 

Current  value  of  missile  lateral 

displacement 

VARTB  (9) 

Current  value  of  missile  altitude 

VARTB  (10) 

Current  value  of  missile  pitch  angle 

VARTB  (11) 

Current  value  of  missile  yaw  angle 

VARTB  (12) 

Current  value  of  missile  roll  rate 

VARTB  (13) 

Current  value  of  missile  roll  angle 

4 
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TABLE  3.  PROGRAM  CONSTANT  INPUTS 


Fortran  Name 

Definition 

Units 

THETO 

Initial  pitch  angle 

rad 

PS  10 

Initial  yaw  angle 

rad 

PHIIN 

Initial  roll  angle 

rad 

WTIN 

Initial  missile  total  weight 

lb 

TT 

Current  value  of  time 

sec 

WO 

Initial  velocity  component  of  missile 
in  pitch  plane 

ft/ sec 

U 

Initial  velocity  component  of  missile 
along  the  longitudinal  axis 

ft/ sec 

VO 

Initial  velocity  component  of  missile 
in  yaw  plane 

ft/ sec 

RO 

Initial  yaw  rate  of  missile 

rad/ sec 

QO 

Initial  pitch  rate  of  missile 

rad/ sec 

PROLL 

Initial  roll  rate  of  missile 

rad/ sec 

XR 

Initial  range  of  missile  in  earth 
coordinate  system 

ft 

Y 

Initial  lateral  displacement  in 
earth  coordinate  system 

ft 

Z 

Initial  vertical  displacement  in 
earth  coordinate  system 

ft 

XX 

Number  of  time-steps  per  print  out 

- 

XXTIN 

Time  of  first  print  out 

sec 

TIMP 

Total  impulse 

lb- sec 

TPROP 

Total  propellant  weight 

lb 

DELT 

Time  integration  step 

sec 

XII 

Initial  transverse  moment  of  inertia 

slug- ft2 

XIB 

Final  transverse  moment  of  inertia 

slug- ft2 

TTO 

Initial  time  for  table  values 

sec 

TLI 

Initial  thrust  lever  arm 

cal 

TLF 

Final  thrust  lever  arm 

cal 

WY 

Boost  phase  lateral  wind  component 

ft/sec 

WYl 

Coast  phase  lateral  wind  component 

ft/ sec 
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TABLE  3 . (Continued) 


Fortran  Name 

Definition 

Units 

WX 

Boost  phase  head  wind  component 

ft/ sec 

WX1 

Coast  phase  head  wind  component 

ft/ sec 

TFLD 

Time  of  fin  folding 

sec 

DELMIN 

Minimum  time  integration  step 

sec 

ROLLA 

Initial  value  of  roll  moment  of  inertia 

2 

slug- ft 

ROLLB 

Final  value  of  roll  moment  of  inertia 

2 

slug-ft 

CLAI 

Initial  value  of  control  lever  arm 

cal 

TB 

Boost  burn  time 

sec 

I BO 

Index  indicating  whether  to  print 

0 - no  print 

out  angular  errors  after  burn 

1 - extra 
print 

BL 

Missile  body  length 

ft 

TMAX 

Maximum  time  - simulation  stops  after 
this  time 

sec 

LCS 

Index  indicating  last  case  of  dac/i 

O-not  last 
case; 

1-last  case 

QE 

Quadrant  elevation 

rad 

S 

Missile  reference  area 

ft2 

B ■ 1 

Thrust  misalignment  angle 

rad 

PHIO 

Initial  thrust  misalignment  direction 

rad 

DX 

Missile  diameter 

ft 

ELL 

Effective  launcher  length 

ft 

DCANT 

Fin  cant  angle 

rad 

PCD 

Percent  effective  drag  coefficient 
utilized  in  error  budgets 

CNAI 

Constant  value  of  CNA  before  fins 
are  unfolded 

/rad 

SLI 

Constant  value  of  SL  before  fins 
are  unfolded 

cal 
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TABLE  3.  (Continued) 


Fortran  Name 

Definition 

Units 

T 

Time  table,  independent  variable  for 
thrust  table 

sec 

A 

Thrust  table,  f(T) 

lb 

TTA 

Time  table,  independent  variable  for 
PSPIN  table 

sec 

AR 

PSPIN  table,  f (TTA) 

rad/ sec2 

TTA 

Time  table,  independant  variable  for  ASPIN 
table  (Both  ASPIN  and  PSPIN  will  not  be  used 
in  same  simulation.) 

sec 

AK 

ASPIN  table,  f (TTA) 

rad/ sec 

UA 

Mach  number  table  independent  variable 
for  SL  table,  power  off 

- 

AJ 

SL  table,  power  off,  f(UA) 

cal 

UD 

Mach  number  table,  independent 
variable  for  CD  table,  power  off 

— 

AM 

CD  table,  power  off,  f(UD) 

— 

UJ 

Mach  number  table,  independent  variable 
for  CNA  table,  power  off 

— 

AT 

CNA  table,  power  off,  f(UJ) 

/rad 

UC 

Mach  number  table,  independent  variable 
for  CD  table,  power  on 

— 

AH 

CD  table,  power  on,  f(UC) 

— 

UI 

Mach  number  table,  independent  variable 
for  CLO  table 

— 

AS 

CLO  table,  f(UI) 

— 

UG 

Mach  number  table,  independent  variable 
for  CLP  table 

— 

AO 

CLP  table,  f(UG) 

(rad 

UB 

Mach  number  table.  Independent  variable 
for  CNA  table 

AG 

CNA  table,  f(UB),  power  on 

(rad 

UE 

Mach  number  table.  Independent  variable 
for  CMQ  table 

— 




Fortran  Nane 


THETS 
PS  ITS 


THRUST 


TABLE  4.  PROGRAM  OUTPUT 


Definition 


Current  value  of  time 

Component  of  missile  velocity  along 
longitudinal  axis 

Component  of  missile  velocity  in 
yaw  plane 

Component  of  missile  velocity  in 
pitch  plane 

Pitch  angle 

Yaw  angle 

Roll  angle 

Missile  pitch  rate 

Missile  yaw  rate 

Missile  acceleration 

Velocity  along  the  z-axis  in  the 
earth  coordinate  system 

Velocity  along  the  y-axis  in 
earth  coordinate  system 

Steady-state  pitch  angle  after  burn  out 

Steady-state  yaw  angle  after  bum  out 

Range  of  missile  in  earth  coordinate  system 

Control  force  for  system  with  directional 
control 

Missile  current  thrust 

Lateral  wind  velocity  component 

Pitch  angular  error  at  burn  out 

Yaw  angular  error  at  bum  out 

Velocity  along  the  x-axis  in  the  earth 
coordinate  system 

Missile  flight  path  angle 

Missile  angle  of  attack  in  the  vertical 
plane 


ft/ sec 


ft/ sec 


ft/ sec 


rad/ sec 

rad/ sec 
2 

ft/ sec 
ft/ sec 


ft/ sec 


ft/ sec 


ft/ sec 


E.  Flow  Chart 


SUBROUTINE  DESUB 


III.  PROGRAM  LISTING 


i 

I 

I 

1 


PROGRAM  HAIMINPUT.OUTPUr.TAPf S= INPUT .TAPt6=0UTPUT> 

DIMENSION  VaPTp ( i j) .««« ( J) 

DIMENSION  T (A)  .»  (4)  .UA(13»  . « I < 1 3 ) « UP ( ] 3 > • AG ( 1 3 ) .0C(tJ>.MM(h>. 

**  J (13)  .UG  (2  ) .A</(2)  .Uh  (2)  .AP(2)  .TT A (4)  .A«  (4»  «AK  (4)  «I/T  (7)  .AS  (7*  . 

•0  J (13)  .AT  (13)  «ZAL  (A)  .AUU)  «1Tn(2>  tAL  (?)  .UC  (A ) .flM(*»l  .lit  (?»  .AN  (2) 
COMMON  GWA V »Ut  *P.  3 »CuA  . AN.  « T .0  T v *PH  1 *PP  IO.XK .CD.ttV. SL • X I .DX.CMU 
Common  TL.CL4..-0LL)  .aPh!  .»•  LU  .NELL .MU. MU.P5PIN  .CLP.CLD.DCAN'T »PDOT 
common  Pcn.ASPiA .glo 

tXTEkNAL  iltSLlM 

NAME  LIST /NaMA/1 pETo.PPI iW./.T.TI.U.XK.aT IN .00 • MO .NO. VO .PS 10* 

1 xXTIN.nTM.MY.t-LL.'JE.S.HHl  a.uX.TLI.TLF.POLLA.KX.TP.PPOLL 

<?XI I .Xlb.TlMP.TMROP.TTO.CLAI .CELT .TMAX .T.a.UA, A I ,UH« AG.CC. Ah. A J» 

3TT  A » AL  . UP.  AM  .lit  . *N  «DtLM i u.^OLLH  « IPO.BX  .LCS.QCanT.UG.  AO  .(JH«  AP« 

4TTA.AM.MYl.MXl.PCn.AK,Ul.AS*UJ.AT.TEC0.7AL*AU*HL.CNAl.SLl 

PLOU*bL/l)X 

LCS=0 

C ATMOSPHERIC  DAU  AT  SfcAUtVFL 

TfMP0=Sltt.6Ptt 
°PSS0=21 16.217 
GMAV0*32. 17404 
BF»GPAVO«. 1631653 
100  MEAD (S.NAMaI 
MMITf  (6.10<i7) 

1097  FOMMAT(lHl) 

•»M1TE(6.NAMA) 

C INITIAL  CONDITIONS 

TpFTA=UE«57.S 
1C*0 

OA»TH(l)=TT 
VAPTB (?) *M0 
VAPTB(l)eU 
V»WT6(*)=V0 
VANTH(S) sPO 
VAR1B(6)>('0 
VuRTR (7) *Xh 

V nWTb I A) xY 

VABT8 (9) *7 

V49TiAU0)=TMs.TV 

VmRTH111)*PSI« 

VABTR(1?»*mR0LL 

VAf»T«(13)*PMll:, 

•TXiTIN 

XXT*XXTIN 

«.LL=1 

C atmospheric  OATa  CALCU-MTfcr. 

4 i'.aD=1.-VAPTM(9)/2«'p5S531.S 
Gh  aV*VjAAVO*  ( I ./DA)4)  •*? 

ALT«-VAPTM(9)/I)Ai) 

U (ALT-36.  »9.2?am>  AS. 45. ah 
AS  TFMP»TuMPO*  (“»li|)3R661t*ALT ) 

PhCS*MHSSO«(ltMD/TtMP0)**>“- 

P*2a.46N0*PPS«:/(<.O714.6/*S7*TEMP) 

Om=  (Temp*?4,0J.07S)«*.S 
•;0  TO  ">0 

4A  (r  M^slM<..9ri7M 

m^ss=4 /?.H765«r »-(GPtVO«(-.UOl)OOl49)*(ALT-360P9.?3R8) ) 
m*?m.RoA  j«pRS5/I9?40i A7.44 
Ow  sSA»H#r7^ 

T?u  * V = *Y 

l)  (TT .uT . Ip ) «v=wY 1 
c A.pOimANIC  04  f 4 CALCDLATrD 

I)  (TT.GT.TR)«X=A*( 

*O*t»X*CCS(VAPT.-l()0)*UF> 
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- 


i 


. 


•bAbXASINTVAMiMlOMu)  ) 

VT0T*S<1RT  I TVARTP  » >/*wf  *i  <-)  -»\i  I '♦  I V»iO|h  ( ?) ) *»2> 

ALP*  T VAPTH  4 2 > -** ) / /TO  T 

ALPHA* SORT  ( TVAWTH<2)-*A),,'»  TyAWIH  !<.»-■  v>**?>/VTuT 

UM-VTOT/OM 

ALTFT*-VARTH19> 

CALL  AIM  TT.A.TT.xT) 

CALL  X1KHTTA.A-.TT.PSPI  .) 
call  xiktitta.aa.tt.wspim 

CALL  XIMTUG.AO.UM.CLP) 

call  xintiuh.ap.iik.clc) 

CALL  XINTTTTb.AL.TT.XA) 

CALL  X1NT  TOC.AK.ltM.CPO) 

Call  XINT(UH.AG,U«.CM«) 

CALL  XlNT  TUI .AS.UM.CLO) 
call  xint tzal.au. ALTFT.AKL) 

RN*RNLWUWHL 

CDf*.ZRA»MLOt)AT).03<»2-.10Sl*UP».J03S*'*P*,?)/TPN**.?) 

IKTTT.GT.Tb)  r.O  TO  lol 
C PHYSICAL  CHARACTERISTICS  UoCATfcC 
G SP1NR«TINP/TPR0P 
■T00T»XT/SPINP 
•T*NT-*TDOT«OCLT 
bBR*4XlB-XII»/TTH-TT0) 

Xi*K J I aBBPaTT 

fH»T*  TTLF-TLI ) / TT0-TTO) 

TL*TU«B»»T*TT 
CLAaAbS TCLAI-TL) 

ROLLI*NOLLA. T tTT-TTO) / T TH-tIO> )• TROLLR-ROLL A) 

CALL  XlNT  IUA.AI.IM.SL) 

CALL  XlNT  TUC.AH.IIN.CU) 

GO  TO  S 

101  CALL  XlNT  (OA.Aj.OM.SL) 

CaLL  XlNT  TUU.AP.ttM.CD 
CALL  XINTTUJ.AT.UP.CKA) 

S XP>fcT/GRAV 
CO*CO*COK 

IFIIT.bUTFLO)faO  TO  ? 

CNA*CnAT 
SL*SL1 
T CONTINUE 

C IKTfGPATION  SUHROOTlNf 

K TGAMTN  T 1 3)  .fit  .0.2832 ) V ART*- ( 1 3 ) « VAPTH  ( 1 3 ) -G*  2H32 
lFT»AR|8T13).Lr.-G.2«32) YArTBtI3)=VAOTSTl3).6.?832 
CALL  SSIHOCTWAPTB.bbM.OCLT .CCLPIN.  1 2. TC.UtSL'B) 

TTaVAPTG T I ) 

AXmVAPTR  T3)*C0S  T YAHTH  (10)  .ft ) ♦ VAPTH  T *)  *5 IM  VAOTfl  ( 1 0 ) *Ofc ) * 

•CoS TVAPT8 1 1 1 ) ) -V API b T* > «S 1 » T V AmTh  T 1 1 > > 

ATaYAHTR TA) .VAMTb ( 3) *SIM  VaM*  T 1 1 > > 

A?*YARTRT2)*C0S  (VARTHTlu)  *<;t  )-VA-T^(l)«SI''  (VAPTH  < 1 0 » *<5C ) 
TH£TA*4TAM T-AZ*/AX)«57.? 

T)t.AG*.S«R*S*Cf)»T  VAPTH  T3)«*2) 

ACC*  TXT—DhAG)  /XP 
If TTT.LT.TH)GO  To  200 
C STEADY  STATE  ANGLES 

THETS*VARTH(JA»-V!\PTH(2)/VAPlHT3)-WAOTrT6)**I/Tnx*XP«VAPTHn)t*SL). 
1 TVAKTfi!2)-b?)»P*Sn>X*CN'5/U.«AP«VARTH  C*)*SL) 

PSITS* VARTh ( IT ) «VAPTfc (A)/(vAKlPT3>«C0S(Rt).VAfiTHT2)#SIKT')K))- 
IYaRTNT5)*XI/TPX«XM*VaPThT3»«SL*C0S<0E)  ) - ( VaPTh  ( A ) -pY  ) *P*S*i)X*C.TU/ 
2TA.axm*VARTFi(3)*SL«C0ST0E)  > 

IF  TTT.tof .XXTIPMITP ( A. 1 OR®) THtIS.PS ITS 
200  IKtlnO.EQ.OGO  TO  201 
C AK6ULAW  tPPOHS  AT  BOOST  HUrK.  OUT 

V2*VAPTHf?>*C0Sl0AM«(lo>  )-VAPTH)3)*STMVAPTB(10) ) 

YT*YA*TH  T3)  *COS  ( VAKTo  ( 1 .1)  ) »S  It>  (WaRT»(1  I ) ) »vAHTh  )«. ) *COS  T VAPTH  ( 1 1 ) ) 
l«VAPTHT2)*SINTVAPTbT10) )*STK (VAuTriTll) ) 
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► 


E 


ZfeW**V2/VT0T 

rtB6*vr/vToT 

IF  <TT.b£.XAT)W4ITt  (6.104H)  yZ. vr.ZEPB.YE4* 

2ul  IF (TT-XXT) Vl .90.90 

90  6hIT€(6.109G)VAi;T4<l).VABT4<3).VAHTB<7).AP 

4BITE  (6*  1492  > VAST« ( 1 1 ) . VA^TH  (4)  .AY.VAPTH  <H) 
»*ITt"(6.1(;93)VA«;TP<13)»V»kTr<6).VABTH<S)  .UBA6 
««IT£  <6.1C9‘>)XT.'.r.XS,«C< 

«MTE(6.11O0)  AX.Th-Ta.ALP 
U*IT? (6.1096) 

109U0F0BPAT (8H  TIME*  *tl6.7*2«.«HOs  «€  16.  7«2X.4*iAB*  .£14.7.2*. 

17FPASS*  .E14.7) 

10910FObp*T(8H  TmET*  «E 14.7.2A.4~«x  .£14. 7.?x.4haz=  .£14. 7. 2*. 

17h?»  .tl4.7) 

1092flFv>BHAT(8M  PSIO»  .£14.  7,2*.4HVx  .t  l6.7,2da6H»V=  .€16.7.2*. 

17MY*  .EI4.7) 

10930F OBPAT (PM  PMI*  .L14.7.2A.4HC*  »£  14.  7,2X  .4PB=  .£14.7.2*. 

l/H)BAGa  .€14.7) 

1095  FOOPAT  (8H  T»iMIST«.tl4.7.2A«4na7*  ♦ £ |4. 7.21  .a*»Xk  = .£14.7«2X.4HACC*. 
If 14.7) 

1096  FOBPATllM  ) 

109®  FoPPATlAM  V2«  •£14.7.2X.4HVY*  .€  14„7*2X»4Fi7£x  .E14.7.2X.4HYE*  . 

If  14.7) 

1099  F0RPAT18M  THE7S*  .E14.7.2X,4hPTS«»£14.7) 

1100  F©«PAT(8H  A007=  .l14.7.?a«4hThAx,ei4.7.2X.4HAlP*.E14.7) 

XXT>XXT»XX*OELT 

91  CONTINUE 

C PN06MAM  STOPS  04  NEGATIVE  ALTITUDE  ©«  PAX  TfPF 
IF (VANTH (9) .6T. 1 0. >G0  TO  30v 
IF(TT-TPAX)4.30r..7u0 

300  «MTE<6«lV<0IVA4Tt<<l>  .VAt»T8(7>  .*P 

•BITE <6»1P91 1 VABTB (IS) « VAmTE (?) .A2.VAPT8 (9) 

•FITE  <6.10921 VABTn  (II).  VA4TI!  (4)  .AY. VA9TH  <*> 

•**ITE  (6.1.193)  VA4  (4(131  «VA*Tt!(6>*VAWTH(S)  »C»AG 
>*6|TF  (6.109b)  Xl.aV.Xiv.AU 
661T€(6.109A)V2.VY*2fc»k.YEc4 

C LAST  SET  OF  DATA  MGNALLtu  teY  LCS-I 

301  IF (LCS.EO.O)GO  TO  IOC 
C«LL  EXIT 

t\0 

SUP40DTINE  DFSOP(V.F.J) 

COPPON  GH4V.WE .P.S.CNA »*)•»» 1 .DTP .PHI .PHlO.XK.CD.OV.SL.XIaDH.CMO 
COPPOFI  TL.CLA.kOLLi.APHl.tLL .KELL..P.PD.PSPI6.CEP.CLD.0CANT.P00T 
COPPON  PCO.ASPlx.CLO 
OlPENSION  V ( 13) 

1F(V(7).GT.FLL)«FLL*? 

Go  TO  <10.20. 30.40. So.64.7c.OU.9u.100. I Iv. 120). J 

10  F.o 

GO  TO  <12.11 l.KFLL 

1 1 FaGBA V*COS ( V (10) .PF ) «V  <6)*V (3)— (V ( 3)  *»? ) ap*S«CNA*  ( V (21— (HI)  / 

1 < ?.•*"• (V (3) ~PU) ) • 

2(XT4(>TP«SIN(V(l?)«PH10MAK*V(lt)  )/XP 

12  kr  TlikN 

20  f s ( XT-.S*P*S*PCr*C')*V (3) **?> /XP.Y(4)*V(S)-V(2)*V(6) 

GO  TO  (22.21).  «tLL 

21  F*(AT-.S*H»S«PLr,Cl).V(3)**?)/AP-GP»v.SI4(V(10).QE)*V(4)*V(S)-V(2)4 


1 V < 6 ) 

22  P*  TUP:. 

JO  F*0 

04  TO  (32. Jl). SELL 

31  F * (“l.*V(S)*V<3))-(V (4)-.V>* (V 13) *•?) •P*S*CNA/(?.*XP4< V(3)— WU) ) • 
) <XT*r>T»*CCS(V(13).kHl<>).Aiv«Y<il ) )/XP 

32  KtTlIFw 
4«  F >U 

ou  lu  (42.4)). KELL 
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41  F«  < V C3 » ••2 » *»€••■!**  (V  VI  I- 

1 .S* (DX*»2) *P*S* V ( S) *V O ) * 

2 C**Q/<2.»X1I*<-1  .*xT*CT-‘*C<s£  ( V ( 1 3>  *P«I  o>  ML*D*-»  .*nx*CLA*XK* 
3V(ll))/XI*R0LLi«V(l2>«V  (r)/M 

42  RETURN 

50  F«0 

GO  TO  (52.51 ) .Ki.LL 

51  Fb  (-1.*(V(J)**2>*p*S*I'X*cN “•  (V  (?»--».  )*SL>  /(?.«XI«(V( 3) -MU)  >- 
1 (.5«P«S*(DX»b?>*V<»>>* 

2VO>»C*10)/(2.*XT)  *(D*«TL»*T*,jIP*slN<v(n)-»PRlO>-t)A*XK*v(l..)«CLA>/ 
3X1— ROLL  I*V(12)’*v(5)/Al 

52  RETURN 

(SO  FBV(3>*COS(V(10)*Ct)*V(?)«SlU(VOo>.OE)*COS(V(H>  > -V (A) *SIN ( V < 1 1 > ) 
RETURN 

70  F*0 

GO  TO  (72*71 ) .KELL 

71  F*V(3)»COS(V(10»*O£l»blN(Vm»)*V(?)*51MV(l«)*OE)*51N(V(ll>)* 

IV (4)*C0S ( V (11)) 

72  RETURN 

BO  FbV(2)*COS(V ( 10)  »Qt  >-VI3)*SlN(V(10) ♦V > 

RETURN 

90  F«0 

GO  TO  (92.vlt.KfLL 

91  E*V (6) 

92  RETURN 

100  FbO 

60  TO  (102*101). KFLL 

101  F«4(5)/COS(V(10) *Ot) 

102  RETURN 

110  Fb— R*S*  (OX»*2)  #V  (3)#V(12)  #,'Lr/  (4.  *>»0LL  I ) *F*5*  ( 0 ( 3) ••?) #DX*CLD* 
1DC4NT/(?.*ROLL1)»PSPIN*.S*'»«S«0*»CLO*V(3)»»?/ROU.I 
FbO.  — ■■ 


PUOTbF 
RETURN 
120  F« v ( 12) 
FbASPIN1 
RETURN 
END 
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«m  u mm  oonun  or 


SUBROUTINE 

SOLVCE  (V.m.DEL.PELMIN'.N.IC.OFSUFI) 

SVOEOOOO 

c 

SOLE  00  in 

c 

GIVEN 

AN  N TH  ORDER  DIFFERENTIAL  FOUATION 

SVUE0020 

c 

0<N)r/0X(F>)  b F ( A 

.T.YR.T2R»....Y(N-l)P)  MMERE 

P DENOTES  PRIME  SVDE0030 

c 

SVUEOORO 

c 

1 V 

VARIAaLt 

I^btr  ifet 

SVUEoOSG 

c 

V ( 1 1 

* A 

S VUE 0060 

c 

V (’) 

B V 

SVUEO07C 

c 

VO) 

b Y FftIMF 

SVUt  JOB'S 

c 

V (A) 

= Y T«0  Prime 

SVUE0090 

c 

SVLE0100 

c 

V (N. 1 ) 

* Y (N— I ) Prime 

SVUtOllO 

c 

V(V*?> 

s 1 1 (N  ) V /f)A  (N)  UPON  RFTURN 

SVUt0120 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


? N 


MOCKING  AR»aV.  « DtFFFRtNT  ARRAY  MUST  BE  USED  F OR 
t*CH  DIFFERENT  COUATION  OR  SYSTEM  OF  ECt'ATIONS  BEING 

lvaluateu.  thf  uJhension  rust  re  at  least  «n>. 

STEP  SI/E  FOP  INDEPENDENT  VARIABLE.  AN  INCREASING— 

de Creasing  vapIaple  step  sue  (r>  is  usfu  internally 

IN  T*t  ROUTINE.  RETURN  TO  TH£  CALLING  PROGRAM  MILL 
B£  AT  V ( 1 > *ur * • 

4 OFLMIN  MIM»*IIN  VALUE  To  RE  USED  FOR  (H)  WHEN  DECREASING 
Tr*E  INTEGRATION  INCRtPFNT. 

ORDER  OF  T**€  (fiUATIOF  TO  r*F  SOLVPn. 

(JbFD  F Or  I'vTt-NAL  CONTROL.  TrIS  VARIABLE  RUST  RF 
<r'T  EOUmL  TO  7fR0  FOR  THF  FIRST  INTEGRATION  STEP. 

a r iffe.rent  v ar i ahll  rust  me  used  for  fao-  i>ifff_rEut 

LOCATION  OR  SYSTEM  OF  FOUATIONS  -flNG  EVALUATED.  uO 

not  use  a literal  sincf  ic  is  continually  Changed. 


3 UEL 


S N 
(S  IC 


SVUEOUi 
SVUE014v 
SVUtOlSf 
SVUtolBi 
SVUtOl  7t 

svoeoim 

SVOEOl  V( 

SVUtO?G' 

SVUtO?i 

SVUt022< 

SVU*o23« 

SVUF J?4( 

SVUtO?5( 

SVOE024( 

SVUE027C 


i * 


t 
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OU<JOOOOUUO<JO<JOOUOO<JOOOUOWOUOO<JOOO 


7 UESUb  NAVE  OF  ROOliNt  CONTAINING  IrE  EQUATION  TO  HE  SOLOED. 

T»F  N«Mt  USED  In  T Ht  CALLIN'.*  PROGRAM  MOST  BE  DCF  Inti) 
IN  A EXTERNAL  STATEMENT  IN  IrE  CALLING  PROGRAM. 

SUhhOUT INE  SSIMOE • V.v.OtL • OcLMNtN* IC*OFSU«) 

GIVEN  THE  FIPST  OPOER  SYSTEM  OF  DIFFERENTIAL  EQUATIONS 
DYI/OX  * FI  (X*Y1*...*YM 
DY2/DX  * F»IX.Y1*....YN> 


OVN/DX  r FN  (X.  YI  « . • . • YM 

V VAC  I ABLE  TAMLC  IE* 

VIII  * X . 

V «5»»  * Yl 

VITI  = Vr 


V (N* | ) a YU 

? « SAVE  AS  2 ArOvE. 

3 DEL  SAVE  AS  3 AROvt . 

A DELNIN  SAVE  AS  A -NOvE . 

*>  N NUMBER  OF  SIMULTANEOUS  LOCATIONS  to  me  solved. 

6 IC  SAVE  AS  6 AMOVE. 

7 OESUB  SAVE  AS  7 ABOVE. 

Ttit  VARIABLES  C AND  VS  IN  THIS  PROGRAM  MUST  BE  DIMENSIONED 
AT  LEAST  CU.M  «N0  VS|K«2). 

10-  APRIL  60 
RALPH  SELIEPS 

DOUBLE  PRECISION  V.MtDEL tUFLMlN.C* VS*TST0P.M*H?*F.TMP1 
EXTERNAL  UESUB 

DIMENSION  VII)*V(I>  *C(A.|S)«VS(17) 

1SYS  » 0 
NR l * N«I 
NR2  > NP|«t 
MOT  « NP? 

GO  TO  3 
ENTRY  SSIMDF 
1SYS  a I 
API  a N.l 
NP2  * NPI«1 
MOT  « NPI 

CORPUTt  STOP  FOP  THIS  INTEGRATION  STtP 
3 TSToP  a v<ll  ♦ hel 

CHECK  E OR  FIPST  TIME  I'.To  ROUTINE 
10  IF  I Ic  .Gt.  01  GO  TO  2u 
FIRST  UMr  IN 
1J  H a t)EL 
NR  a 0 
»<?>  * f'El 
GO  T(  GO 

RRSTArT  IF  INPUT  PEL  h»S  CHANGED 
RV  IMPEL  .*«.  V 12)  1 GO  TO  1? 

R a win 
\<-a  i"(J) 

TEST  FOR  OOLBLING  H.  INCREASE.  TEST  AS  NIMPER  Of  HALVIM6S 
r«C«.E**SRS. 

rO  !► < IC  .LT.  lu«(NH*l)>  GO  TO  60 
IL  * 1 

DvcRLE  H 
Sw  RR  * ?.J*R 

1MV(I)*H2  .GT.  T*>  I OR)  GO  TO  60 
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if  1*2  .LE.  DEL ) 00  10 

-3 

SVOE09.C' 

00  TO  60 

SVUt V93y 

S3 

h * M? 

svuto9ou 

c 

INCREASE  COINTEP 

SVUc097i 

*0 

IC  * IC  ♦ 1 

SV0E?9o0 

c 

INTEGRA!  USING  R-K 

SV0E099o 

c 

SAVE  V TARLE  VALUES 

IN  VS  ARRAY 

SVUEIOOO 

10S 

06  106  1*1 .NTOT 

svutui: 

10b 

vsm  * vin 

SVDtl02C> 

c 

FIRST  P*SS  THRU  »-k 

COhRUTaT 10n 

SVOE  1 1130 

110 

00  119  J*1.N 

SVUt  10HIJ 

IF  IISVS  .NE.  0)  GO  TO 

lit 

SVOElObo 

IFIJ  .EO.  N 1 GO  TO 

11' 

SVDE1060 

F * VIJ.2I 

SVOtlO/O 

GO  TO  119 

SVOt 1063 

lib 

CALL  DESUSIVrf »j) 

SVUt 10 VO 

119 

Cll.JI  * F • H 

SVUEUOP 

c 

SECOND  AND  THIRD  PASS  ft- HU  H-K  CA1CULATI0N 

SvOEUlo 

120 

veil  * VSUI  ♦ .S*H 

SVUt 1120 

DO  ISO  1*2*3 

SVOE  1 1 JU 

DO  130  J*2«NP1 

SVOt 11*1 

130 

V<J»  * VSI.JI  « ,5«CII-1.J-1  ) 

SVUt 11  SO 

00  1*0  J*»«N 

SVOEllbO 

If  IISVS  .N£.  01  GO  TO  1JG 

SVUE117o 

IF  4 J .10.  N I GO  TO  13S 

SVOE 11 80 

F * VI J*2) 

SV0E119U 

GU  10  1*0 

SVUt 1200 

13% 

COU.  DESUIMV.r.JI 

SV0E121U 

— 

1*0 

Cll.J)  * F • m 

SVUE1220 

ISO 

CONTINUE 

SV0E1230 

IF|.5*M  ,LT.  PELNINI  GO  TO  160 

S VUE  12*0 

c 

TEST  FOR  HALVING  H 

SVuEl2aC 

00  153  J*1.N 

S VUE 1260 

TPP1  * IC 12. J) -C 13. J> 1 / IC  1 1 » J) “C 12*  Jl  > 

SVOE127A 

It IA6SITMP1)  .GT.  .02SI  GO  TO  1SS 

SVUt 1280 

1S3 

CONTINUE 

SVOt 129u 

£0  TO  160 

SVUE130C 

c- 

HALVE  H 

SVOC1310 

15% 

H * #5*H 

SVOt 1320 

M*  * NM* 1 

SVUE1330 

IC  * 1 

SVUt 13*o 

00  157  J*1.NT0T 

svueijso 

1ST 

VIJI  « VSIJI 

SVOE1360 

GO  TO  110 

SVOE1370 

C 

FOURTH  PASS  THt»U  G-A  CALCULATION 

SvOE 1 3bo 

16u 

VIII  • VSIII  * t 

SVUt 1 390 

RO  161  J«2.NP| 

SVut 1*00 

161 

vui  * vsiji  ♦ cn.j-ii 

SVUt 1*10 

* 

CO  1 70  J*|«N 

SVUt 1*23 

IF IISVS  .NE.  0)  GO  TO  165 

SVut  1 * J(. 

IFIJ  ,L0.  N 1 GO  10  |6s 

SVOE  1**1,. 

/ 

f ■ VIJ.2I 

SVOt  1*3(1 

uO  TO  170 

SVOE 1*60 

Ibb 

CALL  OESUttltf.r.jl 

SVOE 1*?C 

1T0 

C I*.  Jl  « F • M 

SVUE 1*60 

c 

UPDATE  tf  TAHLf 

SVOE 1*90 

00  1*6  J*I «K 

SVUE1S00 

1 8v 

v|J.|l*VSIJ*ll*ICIl.JI  * 2..U0* IC 12* J) *C IT. J) • ♦ CI*.JII/  6.000 

SVOE1S10 

c 

TEST  tOP  f>C  Of  1 N T t Gw a 1 ION  STEP 

SVOE 152o 

200 

if  IV 11 1 • H/2.0  .Gl.  TSTOR»  Gf>  TO  220 

SVOE 1530 

GO  TO  *0 

SVUE 15*0 

c 

Er0  Of  PEGUlRfo  INTEGaaT  I ON 

SVOt  1V»J 

220 

VII)  * TSTOP 

SVOE  1560 

IF IISVS  .Nt.  0)  GV  TO  v«9 

SVOt 157u 

Call  UtSUHIV. VU.P2I.I) 

SVOt 1560 

999 

CONTINUE 

StUtlSV. 

30 

tad) 

4(?> 

•M3) 


H 

DtL 

NM 


10 


00000 


HC TURN 
END 

SUBROUTINE  XINT <XX, YY«KtY ) 

DIMENSION  XX ( IS) • YV ( 15) 

I»1 

IF  (X.LT.XX(I*1))G0  TO  10 
I»I*1 

GO  TO  S 

Y>YV(I)*((X-XX(I))/(XX(I*1)>XX(I)) )*(YY(i«l)-VY(I)) 

“t  TURN 
END 

(OOOOOOoOOOOOOOOO 


SVDE1600 
SVUE1610 
SVDE16£0 
S VUE  1630 
S VUE  1640 
SVOE16SO 
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